!************************************************************************************
!Module containing subroutines used in amoeba (from Numerical Recipes, Appendix C1)
!************************************************************************************
MODULE mod_nrutil

USE mod_types

IMPLICIT NONE

INTERFACE arth
 MODULE  PROCEDURE arth_r, arth_d,  arth_i !arth_k
END INTERFACE arth


INTERFACE assert
 MODULE PROCEDURE assert1, assert2, assert3, assert4, assert_v
END INTERFACE assert


INTERFACE assert_eq
 MODULE PROCEDURE assert_eq2,assert_eq3,assert_eq4,assert_eqn
END INTERFACE assert_eq


INTERFACE reallocate
 MODULE PROCEDURE reallocate_rv,reallocate_rm, reallocate_iv, reallocate_hv, &
        reallocate_im 
END INTERFACE


CONTAINS



!**************************************************************************
!(1) Array function returning an arithmetic progression 
!**************************************************************************
!---------------------
!a. Single Precision
!---------------------
FUNCTION arth_r(first, increment, n)

INTEGER, PARAMETER:: NPAR_ARTH=16, NPAR2_ARTH=8
REAL(sp), INTENT(IN):: first, increment
INTEGER, INTENT(IN):: n !I4B?
REAL(sp), DIMENSION(n):: arth_r
INTEGER:: k,k2         !I4B?
REAL(sp):: temp

IF (n>0) arth_r(1)=first
IF (n<= NPAR_ARTH) THEN
 DO k=2,n
  arth_r(k)=arth_r(k-1)+increment
 END DO
ELSE
 DO k=2,NPAR2_ARTH
  arth_r(k)=arth_r(k-1)+increment
 END DO
temp=increment*NPAR2_ARTH
k=NPAR2_ARTH

DO
 IF(k>=n) EXIT
  k2=k+k
  arth_r(k+1:MIN(k2,n))=temp+arth_r(1:MIN(k,n-k))
  temp=temp+temp
  k=k2
END DO

END IF

END FUNCTION arth_r



!----------------------
!b. Double-precision
!----------------------
FUNCTION arth_d(first, increment, n)

INTEGER, PARAMETER:: NPAR_ARTH=16, NPAR2_ARTH=8
REAL(dp), INTENT(IN):: first, increment
INTEGER, INTENT(IN):: n !I4B?
REAL(dp), DIMENSION(n):: arth_d
INTEGER:: k,k2         !I4B?
REAL(dp):: temp

IF (n>0) arth_d(1)=first
IF (n<= NPAR_ARTH) THEN
 DO k=2,n
  arth_d(k)=arth_d(k-1)+increment
 END DO
ELSE
 DO k=2,NPAR2_ARTH
  arth_d(k)=arth_d(k-1)+increment
 END DO
temp=increment*NPAR2_ARTH
k=NPAR2_ARTH

DO
 IF(k>=n) EXIT
  k2=k+k
  arth_d(k+1:MIN(k2,n))=temp+arth_d(1:MIN(k,n-k))
  temp=temp+temp
  k=k2
END DO

END IF

END FUNCTION arth_d



!------------
!c. Integer
!------------
FUNCTION arth_i(first, increment, n)

INTEGER, PARAMETER:: NPAR_ARTH=16, NPAR2_ARTH=8
INTEGER, INTENT(IN):: first, increment, n !I4B?
INTEGER, DIMENSION(n):: arth_i
INTEGER:: k,k2, temp        !I4B?


IF (n>0) arth_i(1)=first
IF (n<= NPAR_ARTH) THEN
 DO k=2,n
  arth_i(k)=arth_i(k-1)+increment
 END DO
ELSE
 DO k=2,NPAR2_ARTH
  arth_i(k)=arth_i(k-1)+increment
 END DO
temp=increment*NPAR2_ARTH
k=NPAR2_ARTH

DO
 IF(k>=n) EXIT
  k2=k+k
  arth_i(k+1:MIN(k2,n))=temp+arth_i(1:MIN(k,n-k))
  temp=temp+temp
  k=k2
END DO

END IF

END FUNCTION arth_i


!**************************************************************************
!(2) Report and die if any logical is false (used for arg range checking)
!**************************************************************************
!----------------
!a. One argument
!----------------
SUBROUTINE assert1(n1,string)

CHARACTER(LEN=*), INTENT(IN):: string
LOGICAL, INTENT(IN):: n1

IF(.not. n1) THEN
 WRITE(3000,*) "NRERROR: an assertion failed with this tag:", string
 PRINT*, "NRERROR: an assertion failed with this tag:", string
 STOP 'Program terminated by assert1'
END IF

END SUBROUTINE assert1


!-------------------
!b. Two arguments
!-------------------
SUBROUTINE assert2(n1,n2,string)

CHARACTER(LEN=*), INTENT(IN):: string
LOGICAL, INTENT(IN):: n1, n2

IF(.not. (n1 .and. n2)) THEN
 WRITE(3000,*) "NRERROR: an assertion failed with this tag:", string
 PRINT*, "NRERROR: an assertion failed with this tag:", string
 STOP 'Program terminated by assert2'
END IF

END SUBROUTINE assert2



!---------------------
!c. Threearguments
!---------------------
SUBROUTINE assert3(n1,n2,n3,string)

CHARACTER(LEN=*), INTENT(IN):: string
LOGICAL, INTENT(IN):: n1, n2, n3

IF(.not. (n1 .and. n2 .and. n3)) THEN
 WRITE(3000,*) "NRERROR: an assertion failed with this tag:", string
 PRINT*, "NRERROR: an assertion failed with this tag:", string
 STOP 'Program terminated by assert3'
END IF

END SUBROUTINE assert3



!---------------------
!d. Four arguments
!---------------------
SUBROUTINE assert4(n1,n2,n3,n4,string)

CHARACTER(LEN=*), INTENT(IN):: string
LOGICAL, INTENT(IN):: n1, n2, n3, n4

IF(.not. (n1 .and. n2 .and. n3 .and. n4)) THEN
 WRITE(3000,*) "NRERROR: an assertion failed with this tag:", string
 PRINT*, "NRERROR: an assertion failed with this tag:", string
 STOP 'Program terminated by assert4'
END IF

END SUBROUTINE assert4



!----------------
!e. N arguments
!----------------
SUBROUTINE assert_v(n,string)

CHARACTER(LEN=*), INTENT(IN):: string
LOGICAL, DIMENSION(:), INTENT(IN):: n

IF(.not. ALL(n)) THEN
 WRITE(3000,*) "NRERROR: an assertion failed with this tag:", string
 PRINT*, "NRERROR: an assertion failed with this tag:", string
 STOP 'Program terminated by assert_v'
END IF

END SUBROUTINE assert_v



!*******************************************************************************************
!(3) To assert the equality of three numbers (report and die if integers are not all equal)
!*******************************************************************************************
!-------------
!a. Two terms
!-------------
FUNCTION assert_eq2(n1,n2,string)

CHARACTER(len=*), INTENT(in):: string
INTEGER, INTENT(IN):: n1,n2
INTEGER:: assert_eq2

IF (n1==n2) THEN
  assert_eq2 = n1
ELSE
 WRITE(3000,*) "NRERROR: an assert_eq failed with this tag:", string
 PRINT*, "NRERROR: an assertion failed with this tag:", string
STOP "Program terminated by assert_eq2"
END IF

END FUNCTION assert_eq2



!----------------
!b. Three terms
!----------------
FUNCTION assert_eq3(n1,n2,n3,string)

CHARACTER(len=*), INTENT(in):: string
INTEGER, INTENT(IN):: n1,n2,n3
INTEGER:: assert_eq3

IF (n1==n2 .and. n2==n3) THEN
  assert_eq3 = n1
ELSE
 WRITE(3000,*) "NRERROR: an assert_eq failed with this tag:", string
 PRINT*, "NRERROR: an assertion failed with this tag:", string
STOP "Program terminated by assert_eq3"
END IF

END FUNCTION assert_eq3


!----------------
!c. Four terms
!----------------
FUNCTION assert_eq4(n1,n2,n3,n4,string)

CHARACTER(len=*), INTENT(in):: string
INTEGER, INTENT(IN):: n1,n2,n3,n4
INTEGER:: assert_eq4

IF (n1==n2 .and. n2==n3 .and. n3==n4) THEN
  assert_eq4 = n1
ELSE
 WRITE(3000,*) "NRERROR: an assert_eq failed with this tag:", string
 PRINT*, "NRERROR: an assert_eq failed with this tag:", string
STOP "Program terminated by assert_eq4"
END IF

END FUNCTION assert_eq4


!-------------
!d. N terms
!-------------
FUNCTION assert_eqn(nn, string)

CHARACTER(len=*), INTENT(in):: string
INTEGER, DIMENSION(:), INTENT(IN):: nn
INTEGER:: assert_eqn


IF (ALL(nn(2:) == nn(1))) THEN
  assert_eqn = nn(1)
ELSE
 WRITE(3000,*) "NRERROR: an assertion failed with this tag:", string
 PRINT*, "NRERROR: an assert_eq failed with this tag:", string
STOP "Program terminated by assert_eqn"
END IF

END FUNCTION assert_eqn



!**************************************************************************
!(4) To find index of minloc on an array
!**************************************************************************
FUNCTION iminloc(arr)

REAL(dp), DIMENSION(:), INTENT(IN):: arr
INTEGER, DIMENSION(1):: imin
INTEGER:: iminloc

imin = MINLOC(arr(:))
iminloc = imin(1)

END FUNCTION iminloc


!**************************************************************************
!(5) To find index of maxloc on an array
!**************************************************************************
FUNCTION imaxloc(arr)

REAL(dp), DIMENSION(:), INTENT(IN):: arr
INTEGER, DIMENSION(1):: imax
INTEGER:: imaxloc

imax = MAXLOC(arr(:))
imaxloc = imax(1)

END FUNCTION imaxloc



!**************************************************************************
!(6) To exchange the contents of two vectors or scalars
!**************************************************************************
!-------------
!a. Vectors
!-------------
SUBROUTINE swap(a,b)

REAL(dp), DIMENSION(:), INTENT(INOUT):: a,b
REAL(dp), DIMENSION(SIZE(a)):: dum
  dum = a
  a = b
  b = dum

END SUBROUTINE swap


!-------------
!b. Scalars
!-------------
SUBROUTINE	swap1(a,b)

REAL(dp), INTENT(INOUT):: a,b
REAL(dp):: dum
  dum = a
  a = b
  b = dum

END SUBROUTINE swap1


!-------------
!c. Integers
!-------------
SUBROUTINE	swap_i(a,b)

INTEGER, INTENT(INOUT):: a,b
INTEGER:: dum
  dum = a
  a = b
  b = dum

END SUBROUTINE swap_i



!**************************************************************************
!(7) To report a message, then die
!**************************************************************************
SUBROUTINE nrerror(string)

CHARACTER(len=*), INTENT(IN):: string
 WRITE(3000,*) "NRERROR: ", string
 PRINT*, "NRERROR: ", string
STOP "Program terminated by nrerror"

END SUBROUTINE nrerror



!**************************************************************************
!(8) To sort an array 'arr'(sp) into ascending numerical order by Shell's 
!    method (diminishing increment sort): 'arr' is replaced on output by 
!    its sorted rearrangement
!**************************************************************************
SUBROUTINE sort_shell(arr)

USE mod_types

IMPLICIT NONE

REAL(sp), DIMENSION(:), INTENT(INOUT):: arr
INTEGER:: i,j,inc,n
REAL(sp):: v

i = 0
j = 0
inc = 0
n   = 0


n   = SIZE(arr)
inc = 1

DO 
 inc = 3*inc+1
 IF(inc > n)EXIT
END DO

DO 
 inc=inc/3
 
 DO i = inc+1,n
    v = arr(i)
	j = i
	DO
	 IF(arr(j-inc) <= v) EXIT
	  arr(j) = arr(j-inc)
	  j = j - inc
	 IF (j <= inc) EXIT
    END DO
	
	arr(j) = v
	
  END DO
	
  IF(inc<=1) EXIT

END DO	 



END SUBROUTINE sort_shell


!*****************************************************************************
!(9) To sort an array 'arr' in ascending order using quicksort, while making 
!    the corresponding rearrangement of the same-size array slave. The 
!    sorting and rearrangement are performed by means of an index array
!*****************************************************************************
SUBROUTINE indexx(arr,jndex)

USE mod_types

IMPLICIT NONE

REAL(dp), DIMENSION(:), INTENT(IN):: arr
INTEGER, DIMENSION(:), INTENT(OUT):: jndex
INTEGER, PARAMETER:: NN=15, NSTACK= 50
REAL(dp):: a
INTEGER:: n, k, i, j, indext, jstack, l, r
INTEGER, DIMENSION(NSTACK):: istack

n = assert_eq(SIZE(jndex), SIZE(arr),'indexx')
jndex = arth(1,1,n)
jstack = 0
l = 1
r = n

DO
IF (r-l < NN) THEN
 DO j = l+1, r
  indext = jndex(j)
  a = arr(indext)
  
  DO i = j-1,l,-1
   IF(arr(jndex(i)) <= a) EXIT
   jndex(i+1) = jndex(i)
  END DO 
  
  jndex(i+1) = indext
 END DO  
 
 IF(jstack == 0) RETURN
 r = istack(jstack)
 l = istack(jstack - 1)
 jstack = jstack - 2

ELSE
 k = (l+r)/2
 CALL swap_i(jndex(k),jndex(l+1))
 CALL icomp_xchg(jndex(l),   jndex(r))
 CALL icomp_xchg(jndex(l+1), jndex(r))
 CALL icomp_xchg(jndex(l),   jndex(l+1))
 
 i = l + 1
 j = r
 indext = jndex(l+1)
 a = arr(indext)

DO
 DO
  
  i = i+1
  IF (arr(jndex(i)) >= a) EXIT
 END DO
 
 DO

  j = j-1
  IF(arr(jndex(j)) <=a) EXIT
 END DO
 IF (j < i) EXIT
 CALL swap_i(jndex(i),jndex(j))

END DO

jndex(l+1) = jndex(j)
jndex(j) = indext
jstack = jstack + 2

IF(jstack > NSTACK) CALL nrerror ('indexx: NSTACK too small')

 IF(r-i+1 >= j-l) THEN
  istack(jstack) = r
  istack(jstack -1) = i
  r = j-1

 ELSE
  istack(jstack) = j -1
  istack(jstack -1) = l 
  l = i
 END IF

END IF
END DO

CONTAINS


SUBROUTINE icomp_xchg(i,j)

IMPLICIT NONE

INTEGER, INTENT(INOUT):: i, j
INTEGER:: swp

IF (arr(j) < arr(i)) THEN
 swp = i
 i = j
 j = swp
END IF

END SUBROUTINE icomp_xchg

END SUBROUTINE indexx



SUBROUTINE sort2(arr,slave)

IMPLICIT NONE

REAL(dp), DIMENSION(:), INTENT(INOUT):: arr, slave
INTEGER:: ndum
INTEGER, DIMENSION(SIZE(arr)):: index
ndum = assert_eq(SIZE(arr), SIZE(slave), 'sort2')

!make the index array
CALL indexx (arr,index)
!sort arr
arr = arr(index)
!rearrange slave
slave=slave(index)

END SUBROUTINE sort2



	SUBROUTINE array_copy_d(src,dest,n_copied,n_not_copied)

	USE mod_types

	IMPLICIT NONE

	REAL(dp), DIMENSION(:), INTENT(IN) :: src
	REAL(dp), DIMENSION(:), INTENT(OUT) :: dest
	INTEGER, INTENT(OUT) :: n_copied, n_not_copied
	n_copied=min(size(src),size(dest))
	n_not_copied=size(src)-n_copied
	dest(1:n_copied)=src(1:n_copied)
	END SUBROUTINE array_copy_d


	SUBROUTINE array_copy_i(src,dest,n_copied,n_not_copied)
	INTEGER, DIMENSION(:), INTENT(IN) :: src
	INTEGER, DIMENSION(:), INTENT(OUT) :: dest
	INTEGER, INTENT(OUT) :: n_copied, n_not_copied
	n_copied=min(size(src),size(dest))
	n_not_copied=size(src)-n_copied
	dest(1:n_copied)=src(1:n_copied)
	END SUBROUTINE array_copy_i




   FUNCTION reallocate_rv(p,n)
 
   USE mod_types
   IMPLICIT NONE

	REAL(dp), DIMENSION(:), POINTER :: p, reallocate_rv
	INTEGER, INTENT(IN) :: n
	INTEGER :: nold,ierr
	allocate(reallocate_rv(n),stat=ierr)!
	if (ierr /= 0) call &
		nrerror('reallocate_rv: problem in attempt to allocate memory')
	if (.not. associated(p)) RETURN
	nold=size(p)
	reallocate_rv(1:min(nold,n))=p(1:min(nold,n))
	deallocate(p)
	END FUNCTION reallocate_rv



	FUNCTION reallocate_iv(p,n)

    USE mod_types
    IMPLICIT NONE

	INTEGER, DIMENSION(:), POINTER :: p, reallocate_iv
	INTEGER, INTENT(IN) :: n
	INTEGER:: nold,ierr
	allocate(reallocate_iv(n),stat=ierr)
	if (ierr /= 0) call &
		nrerror('reallocate_iv: problem in attempt to allocate memory')
	if (.not. associated(p)) RETURN
	nold=size(p)
	reallocate_iv(1:min(nold,n))=p(1:min(nold,n))
	deallocate(p)
	END FUNCTION reallocate_iv



	FUNCTION reallocate_hv(p,n)

	USE mod_types
    IMPLICIT NONE


	CHARACTER(1), DIMENSION(:), POINTER :: p, reallocate_hv
	INTEGER, INTENT(IN) :: n
	INTEGER :: nold,ierr
	allocate(reallocate_hv(n),stat=ierr)
	if (ierr /= 0) call &
		nrerror('reallocate_hv: problem in attempt to allocate memory')
	if (.not. associated(p)) RETURN
	nold=size(p)
	reallocate_hv(1:min(nold,n))=p(1:min(nold,n))
	deallocate(p)
	END FUNCTION reallocate_hv




	FUNCTION reallocate_rm(p,n,m)

    USE mod_types
    IMPLICIT NONE

	REAL(dp), DIMENSION(:,:), POINTER :: p, reallocate_rm
	INTEGER, INTENT(IN) :: n,m
	INTEGER:: nold,mold,ierr


	allocate(reallocate_rm(n,m),stat=ierr)
	if (ierr /= 0) call &
		nrerror('reallocate_rm: problem in attempt to allocate memory')
	if (.not. associated(p)) RETURN
	nold=size(p,1)
	mold=size(p,2)
	reallocate_rm(1:min(nold,n),1:min(mold,m))=&
		p(1:min(nold,n),1:min(mold,m))
	deallocate(p)
	END FUNCTION reallocate_rm



	FUNCTION reallocate_im(p,n,m)
 
    USE mod_types
    IMPLICIT NONE


	INTEGER, DIMENSION(:,:), POINTER :: p, reallocate_im
	INTEGER, INTENT(IN) :: n,m
	INTEGER:: nold,mold,ierr


	allocate(reallocate_im(n,m),stat=ierr)
	if (ierr /= 0) call &
		nrerror('reallocate_im: problem in attempt to allocate memory')
	if (.not. associated(p)) RETURN
	nold=size(p,1)
	mold=size(p,2)
	reallocate_im(1:min(nold,n),1:min(mold,m))=&
		p(1:min(nold,n),1:min(mold,m))
	deallocate(p)
	END FUNCTION reallocate_im



END MODULE mod_nrutil
